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Abstract 

We introduce a Z2 noise for the stochastic estimation of matrix inversion 
and discuss its superiority over other noises including the Gaussian noise. This 
algorithm is applied to the calculation of quark loops in lattice quantum chro- 
modynamics that involves diagonal and off-diagonal traces of the inverse matrix. 
We will point out its usefulness in its applications to estimating determinants, 
eigenvalues, and eigenvectors, as well as its limitations based on the structure 
of the inverse matrix. 



Either formally or as a result of numerical practicality, many physical systems, be 
they classical or quantum mechanical, are boiled down to solving matrix equations. 
As the dimension N x N of the matrix M approaches the limit of solving the matrix 
equation MX = S computationally for a source S with dimension x 1, it becomes 
unfeasible to solve it for an S of dimension N x N using the same algorithm. This is 
so simply because it requires N times more computational time as that of solving for 
a column of S. In physics applications, sometimes one needs to compute quantities 
which amounts to solving for the whole matrix S, e.g. calculating the diagonal and 
off-diagonal traces of M~^. This poses a numerical challenge, sometimes a grand 
one. Especially hard is when the dimension of the matrix M grows fast with the 
physical variables of the problem. For example, when M is represented in the space- 
time coordinates, the dimension N grows as where L is the size of the space-time 
dimension. Quark matrix in lattice guage calculation of quantum chromodynamics 
(QCD) falls in this category. A space-time lattice with merely the size of 16'^ x 24 
gives a quark matrix of the dimension 10^ x 10® including the spin and color degrees 
of freedom. While it is durable to calculate the quark propagator, i.e. M~^(x,0) 
for a point source S at with a reasonably small quark mass (e.g. a fraction of 
the strange quark mass) on today's supercomputers, the quark propagator M~^{x,y) 
from any point to any point is certainly unattainable. For calculations of the 2-point 
functions or 3-point functions with direct insertions, one can get by with the help 
of translational symmetry. But there are cases where one can not rely on such a 
help. These include the calculations of quark loops which are space-time or space 
integrations of the fermion propagators. Examples of interest in QCD include the 
quark condensate and the topological susceptibility with the fermion method 0, 
flavor- singlet meson masses which involve disconnected quark loops in the two-point 
functions, notably the U(l) problem., and the vrA^ a term and the proton spin problem 
which involve quark loop contributions in the three-point functions. 

Instead of waiting for the advent of more powerful hardware, there have been 
several suggestions to solve it approximately with the stochastic approach. One is the 
pseudo-fermion method the other is the stochastic estimation with the Gaussian 
noise 0. We shall pursue the avenue of stochastic estimation with various noises to 
see if one kind of noise is better than the others. 

The idea of introducing noises is hardly new in physics. Historically, it was in- 
troduced to account for the Brownian motion with the Langevin and Fokker-Planck 
equations, to compute the time-dependent correlation functions in statistical mechan- 
ics 1^, for the stochastic formulation of quantum mechanics ^ and quantum field 
theory 0. Stochastic approach to estimating the inverse of an A^ x A^ matrix M 
entails the introduction of an ensemble of L column vectors r] = ri^,...,7]^ (each of 
dimension A^ x 1) with the properties of a white noise, i.e. 



where the stochastic average (■ ■ ■) goes over the ensemble of the noise vectors L; 



{V^) = 0, 
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e.g. {'i]ir]j) = j^J2n=iV?Vj' where r/" is the i-th entry in the noise vector n. The 
expectation value of the matrix element M^^ can be obtained by solving for Xj in 
the matrix equations MX = t] with the L noise vectors rj and then take the ensemble 
average with the j-th entry of r] 

E[M-/] = iv.X,) = Y.Mg'{v,Vk) = Mr/. (2) 

k 

which is the matrix element M^-^ itself. The last step was obtained through eq. (p. 

The pseudofermion method which is based on the Gaussian distribution yields 
identical results as the stochastic estimation. But the stochastic algorithm is not 
limited to the Gaussian noise. Any noise which satisfies eq. (|1]) will work and one 
may work better than another. To see how good a noise is, we study the deviations 
from the orthonormal condition in eq. (|1]) which is strictly true for L ^ oo. For this 
purpose, we define the following 2 errors for measuring the efficiency of a noise. The 
first one is the average absolute value of the off-diagonal element and the second is 
the deviation of the diagonal element from unity 

^HiM-iy- (4) 

i 

For large and independent configurations, we expect the central limit theorem to hold, 
i.e. Ci = = 1, 2, where (jj is the standard deviation. We checked and found this 
relation to hold well for the Gaussian noise, noises with double hump distributions like 
rfe~^ /2 ri'^e~^ and the Z2 microcanonical noise 5{\ 7]^ \ —1) for the dimension 
N = 500. The results of ai and ct2 fitted from the range of L = 10 to 150 are given 
in Table 1. 



Table 1: ai and (J2 as obtained from fitting the Ci and C2 for different noises at 
N = 500 with noise configurations in the range of L = 10 to 150. 



noise 




0^2 


Gaussian 


0.78(1) 


1.49(1) 




0.79(1) 


0.43(1) 




0.80(1) 


0.28(1) 


Z2 


0.77(1) 


0.00(0) 



We see from Table 1 that the off-diagonal deviation ai is about the same for all 
these noises and is close to the asymptotic value of J^/n 0. On the other hand, 
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the diagonal deviation a2 which depends on the deviation of the higher moment; i.e. 
(72 = \/ (^^) ~ 1 decreases as the distribution tends to the bi-nodal form. It is the 
largest for the Gaussian case with an asymptotic value of ^/2 and vanishes for the 
Z2 noise. For this reason, we suspected that the Z2 microcanonical noise may work 
better than the other noises considered here [§. In fact, it has been shown recently 
1^ that the variance of a inverted matrix element due to the stochastic estimation is 
composed of two parts 

Whereas the second part is independent of the kind of noise used, the first part is 
proportional to the square of the diagonal error C2 only. Since Z2, or for that 
matter, has no diagonal error, i.e. C2 = 0, it produces a minimum variance. Other 
noises will have larger variances due to the non- vanishing C2. For comparison of 
different noises, we consider the calculation of chiral condensate which is the trace 
the inverse quark matrix M, i.e. (\E'\E') = TrM~^ /V , for a quenched 16'^ x 24 lattice 
at /? = 6.0 with k = 0.148 for the Wilson action. First we used the conjugate gradient 
program to invert a column of the matrix for a particular gauge configuration and 
find J2ky^i[^kiV /[^ii]'^ — Assuming this ratio is true for all the other columns 
based on translational invariance and extending eq. (|^) to the variance of the trace, 
we find that the standard deviation from the Z2 noise is smaller than that from the 
Gaussian noise by a factor of 1.54. In other words, in order to achieve the same level 
of accuracy , one would need a Gaussian noise configuration 2.4 times larger than 
that of the Z2 noise. Since any noise with C2 7^ will need more statistics than the 
Z2 noise to reach the same accuracy, the Z2 (or Z^v) noise is the optimal choice in 
this sense. 

With the error analysis in eq. (^, we now realize where the stochastic inversion 
algorithm would apply. Since the quark propagator has the generic fall off behavior 
e~™'^~^'/(k ~ y\ + C){a))'^ (a is the lattice spacing) in the space-time separation 
with n = 3 for short distances, quark loops that involve traces near the diagonal, 
i.e \x — y\ ^ a, will have large signals. As long as the far off-diagonal contribution 
to the variance (the second term in eq. (|^)) does not overwhelm the contribution 
from [Mj~^]^ (the first term in eq. (refeq. variance), the square of the matrix element 
of interest, the noise to signal ratio will be of the order ^ for the near-diagonal 

traces from eq. (H). For the trace itself, there is an extra factor of due to 
the translational, gauge and rotational symmetries. This is certainly true for the 
case when the quark mass m is not very small. For m — 0, it remains to be seen 
if the inverse power behavior l/(|x — y\ + 0(a))" is steep enough to curb the off- 
diagonal contribution to give a sufficiently small variance for a reasonable L. But 
when one considers the case when |x — ?/| >> a, the signal drops exponentially while 
the error remains constant. Hence, the noise to signal ratio grows exponentially 
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and the application of the stochastic method is invahdated under this circumstance. 
Therefore, one does not expect the stochastic method to be useful for calculating 
general quark propagators. But it is useful for approximating diagonal and near 
diagonal traces of the inverse matrix if the inverse matrix itself is dominated by the 
diagonal and near diagonal terms. 

To illustrate how the noise works in detail, we employ it to invert the quark 
matrix on a quenched 16^ x 24 lattice at /3 = 6.0 for Wilson fermions {k = .148) for 
a particular gauge configuration. We shall report the results of certain diagonal and 
off-diagonal traces corresponding to the scalar, the pseudoscalar, the vector, and the 
axial-vector currents. The point-split currents are used for the vector and the axial- 
vector cases. The accumulated averages for these currents (summed over the spatial 
points on a time slice) for a gauge configuration are plotted in Fig. 1 against the 
noise configuration L. This shows how they approach the equilibrium up to L = 200. 
Also plotted are the histograms for their distributions. We see that in most cases 
they do tend to stabilize for L = 200. If and when L = 200 is sufficient for the 
purpose of estimating these traces with acceptable errors, this algorithm would save 
the computing time by a factor of 5838 as compared to inverting the full matrix 
(dimension = 1.18 x 10^ in this case) with the brute force approach. 

The Z2 noise has bee employed to calculate quark condensate = TrM~^/V 

and the topological susceptibility with the fermion method [Q, i.e. 
X = ^{Tr{'^^M~^)Tr{'^^M~^)) . The preliminary results §| give reasonable errors 
which are a combination of the errors from the stochastic estimation and the gauge 
configurations. Much like the situation with the glueball masses, we found that 
the two-point functions for the disconnected quark loops are too noisy due to the 
exponential fall off of the two-point function e"^"'^* as a function of t. However, the 
disconnected insertion for the three-point function is different. Since there is only one 
quark loop correlated with the nucleon propagator say, the situation is not as bad as 
the two-point function. Preliminary results on the uNa term and the flavor-singlet 
axial charge g\ are encouraging |TD[ and will be reported elsewhere |TI|]. 

It is worthwhile noting that the stochastic estimation is particularly successful for 
the trace (denotes as Re^^ in Fig. 1). As we remarked before, this is due to the 
translational, color and spin symmetries. As a result, the error is proportional to 
1/a//V where N is the dimension of the matrix. With A^ = 1.18 x 10^ in our case 
and the ratio Y^k^ii^kiT li-^iiT — O-^; we predict the error to signal ratio to be 
VI.5 X 10"'^ from eq. (^ for L = 1. This agrees well with the numerical calculation 
shown in Fig. 1. Given this level of accuracy, it is feasible to apply the stochastic 
method to the calculation of the determinant, the eigenvalues, and the eigenvectors 
of the matrix M which might not be feasible with other algorithms. It is well known 
that for a Hermitian matrix M, the density of states can be written as 

p{\) = ~^\\uiImTrG{X + ie) (6) 
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where G is the inverse matrix of X + ie — M. The determinant of the matrix M can 
be given as 

rfetM = e/^(^)''^^'^^ (7) 

Since the stochastic method is most successful in estimating the trace of the 
inverse matrix M for lattice QCD as we have just remarked and it is known that the 
eigenvalues of M are distributed in a reasonably finite range, it would be worthwhile 
exploring the possibility that this could be an efficient algorithm for calculating the 
determinant. Looking for the poles of the density of states in eq. (^) has been 
frequently used as a way to identify the eigenvalues [l^ and the eigenvectors can be 
obtained from the column of the inverse matrix 

Vi ~ lim Im GikiX + ie) (8) 

£^0 



for any column k [O . 



In conclusion, we have proposed an stochastic algorithm for large matrix inversion 
with the optimal Z2 noise. We show that it is particularly efficient for estimating 
traces and near diagonal traces for matrices whose inverses are dominated by the 
diagonal and near diagonal terms themselves. It is applied to calculate quark loop 
correlations in the vacuum and disconnected quark loop insertions in the three-point 
functions in QCD. Noting that it is most successful in estimating traces, we shall 
explore the feasibility of calculating determinants, eigenvalues, and eigenvectors in 
the future. 
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Figure Caption 

Fig. 1 The first and second columns show the accumulated averages of the real and 
imaginary parts of various current loops as functions of the noise configurations L for 
a gauge configuration. The third and last columns give the corresponding histograms 
for their distributions. 
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